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AN  ACCURATE  CURVED  BOUNDARY  TREATMENT  IN  THE  LATTICE  BOLTZMANN 

METHOD 


RENWEI  MEI*,  LI-SHI  LUCd,  AND  WEI  SHYY* 

Abstract.  The  lattice  Boltzmann  equation  (LBE)  is  an  alternative  kinetic  method  capable  of  solving 
hydrodynamics  for  various  systems.  Major  advantages  of  the  method  are  owing  to  the  fact  that  the  solution 
for  the  particle  distribution  functions  is  explicit,  easy  to  implement,  and  natural  to  parallelize.  Because  the 
method  often  uses  uniform  regular  Cartesian  lattices  in  space,  curved  boundaries  are  often  approximated  by 
a  series  of  stairs  that  leads  to  reduction  in  computational  accuracy.  In  this  work,  a  second-order  accurate 
treatment  of  boundary  condition  in  the  LBE  method  is  developed  for  a  curved  boundary.  The  proposed 
treatment  of  the  curved  boundaries  is  an  improvement  of  a  scheme  due  to  Filippova  and  Hanel.  The  proposed 
treatment  for  curved  boundaries  is  tested  against  several  flow  problems:  2-D  channel  flows  with  constant 
and  oscillating  pressure  gradients  for  which  analytic  solutions  are  known,  flow  due  to  an  impulsively  started 
wall,  lid-driven  square  cavity  flow,  and  uniform  flow  over  a  column  of  circular  cylinders.  The  second-order 
accuracy  is  observed  with  solid  boundary  arbitrarily  placed  between  lattice  nodes.  The  proposed  boundary 
condition  has  well  behaved  stability  characteristics  when  the  relaxation  time  is  close  to  1/2,  the  zero  limit 
of  viscosity.  The  improvement  can  make  a  substantial  contribution  toward  simulating  practical  fluid  flow 
problems  using  the  lattice  Boltzmann  method. 

Key  words,  kinetic  method,  lattice  Boltzmann  equation,  Navier-Stokes  equation,  second  order  boundary 
conditions 

Subject  classification.  Fluid  Mechanics 

1.  Introduction.  There  has  been  a  rapid  progress  in  developing  and  employing  the  method  of  the 
lattice  Boltzmann  equation  (LBE)  [24,  20,  5]  as  an  alternative  computational  technique  for  solving  complex 
fluid  dynamic  problems  (see  the  comprehensive  reviews  in  [3,  6]).  In  a  traditional  method  for  computational 
fluid  dynamics  (CFD),  the  macroscopic  variables,  such  as  velocity  u  and  pressure  p,  are  obtained  by  solving 
the  Navier-Stokes  (NS)  equations  [28,  9,  31].  The  lattice  Boltzmann  equation  approximates  the  kinetic 
equation  for  the  single  particle  (mass)  distribution  function  /(x,  £,  t)  on  the  mesoscopic  level,  such  as  the 
Boltzmann  equation  with  the  single  relaxation  time  approximation  [4]: 

(1.1)  dtf+t-vf  =  -![/-  /(0)], 

where  £  is  the  particle  velocity,  / ^  is  the  equilibrium  distribution  function  (the  Maxwell-Boltzmann  distri¬ 
bution  function),  and  A  is  the  relaxation  time.  The  right  hand  side  (RHS)  of  Eq.  (1.1)  models  the  effect  of 
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Fig.  1.  Two-dimensional  9-velocity  (or  9-bit)  lattice. 


the  fluid  viscosity  on  the  molecular  level  through  the  collision  (relaxation)  process.  The  macroscopic  quanti¬ 
ties  (such  as  mass  density  p  and  momentum  density  pu)  are  the  hydrodynamic  moments  of  the  distribution 


function  /: 

(1.2a) 

P=  j  f(x,  t)d£ , 

(1.2b) 

pu=  [  £/(*,  £>  t)d£ 

It  has  been  shown  that  the  velocity  space  £  can  be  discretized  into  a  finite  set  of  points  {£a}  without  affecting 
the  conservation  laws  [14,  15,  1].  In  the  discretized  velocity  space  the  Boltzmann  equation  (1.1)  becomes 

(1.3)  dtfa  +  ia  •  V/  =  -j[fa  -  ,  (a  =  0,  1,  2,  ....  8  for  2-D) , 

for  the  distribution  function  of  discrete  velocities  fa(x,  t)  =  f(x,  £a,  t).  The  equilibrium  distribution 
function,  faq\  and  the  discrete  velocity  set  {£a}  can  be  derived  explicitly  [14,  15,  1]. 

For  2-D  square  lattice  shown  in  Fig.  1,  we  use  {ea}  to  denote  the  discrete  velocity  set,  and  we  have  [29]: 


(1.4) 


(0,  0), 

(cos[(a  —  l)7r/4],  sin[(a  —  l)7r/4])c, 
(cos[(a  -  1)7t/4],  sin[(a  -  1)7t/4])  \flc , 


a  =  0 , 

ol  =  l,  3,  5,  7, 
a  =  2,  4,  6,  8, 


where  c  =  Sx/Sty  Sx  and  St  are  the  lattice  constant  and  the  time  step  size,  respectively,  and 
(1.5)  fieq}  =wap 


3  9  .  .,  3 

l  +  -^ea-u+^(ea-u)  u 


where 


(1.6) 


a  =  0 , 

Wa  ~  ^  a  =  l,3,  5,  7, 

[  i  a  —  2,  4,  6,  8. 
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With  the  discretized  velocity  space,  the  hydrodynamic  moments  are  given  by: 


(1.7a)  P  =  E/“  =  E^q)’ 

a  a 

(1.7b)  pu  =  E eaU  =  E • 

a  a 

The  speed  of  sound  of  this  model  is  c$  =  c/\/ 3,  and  the  equation  of  state  is  that  of  an  ideal  gas, 


(1.8) 


P  =  C2sP- 
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Equation  (1.3)  is  one  of  numerous  ways  to  model  the  transport  equation  of  /,  Eq.  (1.1). 

Based  on  the  Chapman-Enskog  analysis,  the  solution  for  /a(sc,  i)  may  be  expanded  in  the  form  of 

(1.9)  fa(x,  t)  =  fieq){x ,  t)  +  ef£](x,  *)  +  •••, 

where  e  =  Kn  is  the  formal  expansion  parameter  and  Kn  is  the  Knudsen  number,  which  is  the  ratio  between 
the  mean  free  path  and  macroscopic  flow  characteristic  length.  Substitution  of  Eq.  (1.9)  into  Eq.  (1.3)  with 
a  factor  1/e  in  the  collision  term  leads  to 

(1.10)  /«(*,  0  =  -A  [dtfieq)(x,  t)  +  ea-  Vf^(x,  t)]  . 

Proceeding  with  the  Chapman-Enskog  analysis,  it  can  be  shown  that  the  Euler  equations  can  be  recovered 
from  the  solution  for  faq\  and  the  NS  equations  are  recovered  in  the  near  incompressible  limit  (i.e.,  the 
Mach  number  M  =  \u\/c8  1)  by  the  first  two  terms  in  Eq.  (1.9).  The  viscosity  of  the  fluid  is 

(1.11)  v  =  (?s  A. 

Equation  (1.3)  can  be  further  discretized  in  space  and  time.  The  completely  discretized  form  of  Eq.  (1.3), 
with  the  time  step  St  and  space  step  eaSt ,  is: 

(1.12)  fa{xi  +  ecA)  ^  d"  $t)  “  fa(xii  f)  ~  ~~[fot{xii  ”  fa  0]  > 

where  r  =  A /<S4,  and  x*  is  a  point  in  the  discretized  physical  space.  The  above  equation  is  the  lattice 
Boltzmann  equation  [24,  20,  5]  with  Bhatnagar-Gross-Krook  (BGK)  approximation  [4].  The  left-hand  side 
(LHS)  of  Eq.  (1.12)  is  physically  a  streaming  process  for  particles  while  the  RHS  models  the  collisions  through 
relaxation. 

Although  the  lattice  Boltzmann  equation  historically  originates  from  the  lattice  gas  cellular  automata 
[10,  30],  it  is  indeed  a  special  finite  difference  form  of  the  continuous  Boltzmann  kinetic  equation,  z.e.,  the 
LHS  of  Eq.  (1.1)  is  discretized  along  the  direction  of  the  characteristic  line  with  discretization  of  phase  space 
and  time  tied  together  [14,  15].  The  leading  order  truncation  error  of  such  a  discretization  is  then  taken  into 
account  exactly  by  modifying  the  viscosity  in  the  NS  equation  derived  from  Eq.  (1.12)  to 


The  positivity  of  the  viscosity  thus  requires  that  r  >  1/2.  The  lattice  Boltzmann  scheme  consists  of  two 
computational  steps: 

(1.14a)  Collision:  /a(cCj,  t )  =  fa(xii  f)  —  ~[ fct{xii  f)  ~  fa  ) 

(1.14b)  Streaming:  fa(xi  +  ea8ti  t  +  5t)  =  fa(xi ,  t) , 

where  fa  and  fa  denote  pre-  and  post-collision  state  of  the  distribution  function,  respectively.  The  advantages 
of  solving  the  lattice  Boltzmann  equation  over  the  NS  equations  can  now  be  seen.  In  the  kinetic  equation 
for  fa  given  by  Eq.  (1.3),  the  advection  operator  is  linear  in  the  phase  space  whereas  the  convection  term  is 
nonlinear  in  the  NS  equation.  In  traditional  CFD  methods,  the  pressure  is  typically  obtained  by  solving  the 
Poisson  or  Poisson-like  equation  derived  from  the  incompressible  NS  equations  that  can  be  time  consuming. 
In  the  LBE  method,  the  pressure  is  obtained  through  an  extremely  simple  equation  of  state  p  =  c2sp.  This 
is  an  appealing  feature  of  the  LBE  method.  The  discretized  Eq.  (1.12)  for  fa  is  explicit  in  form,  easy  to 
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Fig.  2.  Layout  of  the  regularly  spaced  lattices  and  curved  wall  boundary .  The  thick  curve  marks  the  boundary  location. 
The  solid  circles  (•)  mark  the  positions  where  particle-boundary  collision  occurs.  The  empty  (o)  and  shaded  (•)  circles  are 
fluid  sites  and  solid  sides,  respectively. 


implement,  and  natural  to  parallelize.  The  collision  step  is  completely  local.  The  streaming  step  takes  very 
little  computational  effort  at  every  time  step. 

However,  unlike  solving  the  NS  equations  for  which  the  non-slip  condition  for  u  on  a  solid  wall  is  satisfied 
at  the  macroscopic  level,  there  is  no  corresponding,  physically  based  boundary  condition  for  fa  on  a  solid 
wall  at  the  mesoscopic  level.  For  a  lattice  node  located  on  the  fluid  side  at  Xf ,  as  illustrated  in  Fig.  2, 
Eq.  (1.14b)  clearly  indicates  a  need  for  the  information  of  at  Xb  on  the  solid  side.  Therefore  all  the  effort  in 
the  previous  treatment  of  the  boundary  conditions  in  the  LBE  models  is  much  focused  on  the  calculation  of 
fa  moving  from  the  wall  into  the  fluid  region.  In  previous  works  of  the  LBE,  the  most  often  used  boundary 
condition  on  the  wall  is  the  so-called  bounce-back  scheme  [32,  11,  19].  In  the  bounce-back  scheme,  after  a 
particle  distribution  fa  streams  from  a  fluid  node  at  xj  to  a  boundary  node  at  x&  along  the  direction  of  ea, 
the  particle  distribution  fa  scatters  back  to  the  node  Xf  along  the  direction  of  (=  -ea)  as  /a.  Since  the 
wall  position  xw  was  forced  to  be  located  at  Xb ,  this  is  referred  to  as  bounce-back  on  the  node  (BBN)  [2]. 
However,  a  finite  slip  velocity  at  the  stationary  wall  exists  [22,  19]  and  the  accuracy  for  the  flow  field  is  thus 
degraded  due  to  the  inaccuracy  of  the  boundary  conditions  [11].  In  simulating  suspension  flows  using  the 
LBE  method,  Ladd  placed  the  solid  walls  in  the  middle  between  the  lattice  nodes  [21].  This  is  referred  to  as 
bounce-back  on  the  link  (BBL).  It  has  been  shown  that  the  BBL  scheme  gives  second-order  accurate  result 
for  straight  walls  [33,  19].  Noble  et  al  developed  a  second-order  accurate  boundary  condition  to  compute 
fa  but  it  is  only  applicable  to  straight  walls  in  triangular  lattice  space  [27].  He  et  al.  generalized  the  scheme 
of  Noble  et  al  to  arbitrary  lattice  [19].  Chen  et  al  placed  the  wall  on  the  lattice  node  so  that  Xb  is  one 
lattice  inside  the  wall  [7].  They  used  an  extrapolation  of  fa  on  the  fluid  side  (including  the  wall  node)  to 
obtain  fa  at  x&.  Zou  and  He  proposed  to  apply  the  BBL  scheme  only  for  the  non-equilibrium  part  of  fa  at 
the  wall  [33]. 

For  a  curved  geometry,  the  use  of  BBL  requires  approximation  of  the  curved  solid  boundary  by  a  series 
of  stair  steps.  The  geometric  integrity  cannot  be  preserved  by  such  an  approximation.  For  high  Reynolds 
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number  flows,  the  integrity  of  geometry  is  important  since  the  vorticity  generation  and  stress  distributions  are 
sensitive  to  the  geometrical  resolution.  To  this  end,  He  and  Luo  proposed  to  use  the  LBE  with  nonuniform 
grid  with  second  order  interpolations  [14, 17, 18].  He  and  Doolen  further  applied  the  interpolation  to  the  LBE 
with  curvilinear  coordinates  or  body-fitted  coordinates  [13].  Mei  and  Shyy  solved  Eq.  (1.3)  in  curvilinear 
coordinates  using  finite  difference  method  [25].  While  the  wall  geometry  is  accurately  preserved  in  body- 
fitted  coordinates,  the  flexibility  to  handle  complex  geometries  is  maintained  by  using  the  numerical  grid 
generation  techniques  common  to  the  Navier-Stokes  solvers.  It  should  be  noted  that  perhaps  the  most 
profound  and  rigorous  theoretical  treatment  of  the  boundary  condition  along  the  wall  is  given  by  Ginzbourg 
and  d’Humieres  [12].  The  scheme  proposed  by  Ginzbourg  and  d’Humieres  is  local  and  accurate  up  to  second 
order  in  Chapman-Enskog  expansion.  However  this  work  has  not  attracted  sufficient  attention  because  its 
implementation  is  not  as  easy  as  the  bounce-back  scheme. 

In  this  work,  a  robust,  second-order  accurate  treatment  for  the  distribution  function  fa  near  a  curved 
boundary  is  developed  based  on  the  method  recently  proposed  by  Filippova  and  Hanel  (hereinafter  referred 
as  FH)  [8].  In  Ref.  [8],  the  boundary  condition  for  fa  on  the  solid  side  is  evaluated  using  Eq.  (1.3)  for  /a, 
and  the  Taylor  series  expansion  in  both  space  and  time  for  fa  near  the  wall.  FH  reported  numerical  results 
for  a  uniform  flow  over  a  cylinder  [8].  However,  it  is  found  in  this  work  that  when  tested  in  a  pressure  driven 
channel  flow  (see  implementation  and  discussions  in  Section  2)  there  is  a  strong  boundary  condition  induced 
instability  when  the  distance  from  the  wall  to  the  first  lattice  on  the  fluid  side  is  less  than  half  of  the  lattice 
size. 

Using  the  Taylor  series  expansion  for  the  velocity  u  near  the  wall,  a  new  treatment  for  fa  near  a  curved 
wall  is  proposed  in  this  work.  While  maintaining  a  second-order  accuracy  of  the  solution  in  handling  curved 
walls,  the  computational  stability  is  improved  so  that  lower  viscosity,  or  higher  Reynolds  number,  can  be 
attained  in  the  LBE  simulations.  The  new  boundary  condition  treatment  is  tested  systematically  to  assess 
the  temporal  and  spatial  accuracy  and  robustness  in  2-D  channel  flow  with  constant  and  oscillating  pressure 
gradients,  flow  due  to  an  impulsively  started  wall,  lid-driven  square  cavity  flow,  and  flow  over  a  column  of 
circular  cylinders.  Detailed  comparisons  for  the  flow  field  are  made  with  either  analytic  solutions  or  well 
resolved  numerical  solutions  of  the  Navier-Stokes  equations  by  using  finite  difference  method.  The  improved 
boundary  treatment  represents  a  significant  step  towards  solving  practically  relevant  flow  problems  using 
the  LBE  method. 


2.  Formulation  for  the  Improved  Boundary  Condition.  Filippova  and  Hanel  considered  a  curved 
boundary  lying  between  the  lattice  nodes  of  spacing  8X,  as  illustrated  in  Fig.  2,  and  briefly  presented  the 
derivation  of  their  scheme  for  the  treatment  of  curved  boundary  in  Ref.  [8].  However,  they  did  not  offer  a 
thorough  explanation  to  justify  the  theoretical  basis  of  their  method.  It  is  instructive  to  first  re-examine 
their  derivation  thoroughly.  Based  on  the  insight  gained,  an  improved  boundary  treatment  is  then  proposed. 


2.1.  Re-examination  of  and  comments  on  Filippova-Hanel’s  treatment.  The  macroscopic  flow 
has  a  characteristic  length  of  L.  The  lattice  nodes  on  the  solid  and  fluid  side  are  denoted  as  x&  and  x/, 
respectively,  in  Fig.  2.  The  filled  small  circles  on  the  boundary,  xw,  denote  the  intersections  of  the  wall  with 
various  lattice  links.  The  boundary  velocity  at  xw ,  the  intersection  with  the  wall  on  the  link  between  x& 
and  X/,  is  uw .  The  fraction  of  the  intersected  link  in  the  fluid  region  is  A,  as  illustrated  in  Fig.  2,  that  is: 

(2.1)  .  A  =  p-ZlA . 

\xf-xb\ 


Obviously,  0  <  A  <  1  and  the  horizontal  or  vertical  distance  between  Xf  and  xw  is  A  •  8X  on  the  square 
lattice.  Suppose  the  particle  momentum  moving  from  Xf  to  x&  is  ea  and  the  reversed  one  from  x&  to  xj  is 
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ea  =  -ea.  After  the  collision  step,  fa  on  the  fluid  side  is  known,  but  not  on  the  solid  side.  (Hereafter  we 
shall  use  ea  and  /5  to  denote  the  velocity  and  the  distribution  function  coming  from  a  solid  node  to  a  fluid 
node,  and  /a  is  the  unknown  to  be  computed.)  To  finish  the  streaming  step, 

(2.2)  fa(x/  =  xb+ea6t,  t  +  6t)  =  fa(Xb,  t) , 

it  is  clear  that  fa(xb,  t)  is  needed.  To  construct  /a(aj&,  t)  based  upon  some  known  information  in  the 
surrounding,  Filippova  and  Hanel  essentially  proposed  to  use  the  following  linear  interpolation  [8]: 

3 

(2.3)  fs(x  b,  t)  =  (  1  -  x)foc(xf,  t )  +  xfaHxb,  t)  -  2wa  p-^(e&-uw) , 

where  uw  =  u(xw ,  t)  is  the  velocity  at  wall,  x  is  the  weighting  factor  (to  be  determined)  that  controls  the 
linear  interpolation  (or  extrapolation)  between  /a(x/ ,  t)  and  fi*\xb,  t),  a  fictitious  equilibrium  distribution 
function  given  by 

T  3  9  3 

(2.4)  /a  0  =  p('C/>  0  1  H"  ^2  ea 'Wft/  +  2^4  (ea ‘^/)  ~  ^  c2Uf  Uf 

In  the  above  equation,  Uf  =  u(x/,  t)  is  the  fluid  velocity  near  the  wall  and  ubf  is  to  be  determined  later.  It 
should  be  emphasized  here  that  the  weighting  factor  x  depends  on  how  ub/  is  chosen.  However,  the  choice 
of  ubf  is  not  unique.  For  example,  either  ubf  =  Uf  or  a  linear  extrapolation  using  ubf  =  [(A  —  l)v.f  +  uw]/A 
appears  reasonable  at  this  stage. 

To  determine  x  in  Eq.  (2.3),  FH  considered  flows  under  the  condition 

(2.5)  J^«l, 

i.e .,  the  flow  has  an  intrinsic  characteristic  time  scale  T  that  is  much  larger  than  the  advection  time  on 
the  lattice  scale,  L/c.  This  “slow-flow”  condition  enabled  FH  to  approximate  /«(*/,  t  +  8t)  in  Eq.  (2.2)  by 
fac{xfi  t)> 

(2.6)  fa(Xf  =Xb  +  ea8t ,  t  +  <**)  =  fa{Xf,  t)  +  St  dtfa  +  ***  - 

For  the  purpose  of  the  order-of-magnitude  estimate,  it  is  seen  that  0(dtfa )  =  0(f&/T)  so  that 

(2.7)  fa(xf,  t  +  St)  =  fa(xf,  t)  1  +  0  (^jj  =fa(Xf,  t)  1  +  O0t^)  «/«(*/»*)• 

It  is  noted  that  under  condition  (2.5)  the  neglected  terms  are  of  0(^*^?)  which  are  much  smaller  than  the 
O(^)  terms  of  the  present  interest  [in  deriving  an  accurate  boundary  condition  for  fa(xb ,  t)\.  Applying 
the  Chapman-Enskog  expansion  in  the  form  given  by  Eqs.  (1.9)  and  (1.10)  and  invoking  the  “slow  flow” 
approximation  defined  by  Eq.  (2.5),  one  obtains 

/«(*/,  t)  =  ft\xf,  t)-\  [dt/<eq)  +  ea  •  V/<eq)]  +  •  •  • 

(2.8)  *  ft\xf,  t)  -  Aea  •  V/ieq)  +  •  •  •  . 

For  given  by  Eq.  (1.5),  the  leading  order  term  in  V/aeq^  is  given  by  3iuapV(iZ'ea)/c 2  since  the  rest  are 
higher  order  terms  in  the  near  incompressible  flow  limit.  Noticing  that  A  =  r^,  Eq.  (2.8)  becomes 

fa(xf ,  t)  «  fieq\xf,  t)  -  3waprSt^e&*V (uf  * e5 ) 

(2.9)  =  /ieq)(®/,  t)  -  6  wap^uf-ea  -  3wapT5t~ea-V(urea) , 
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which  approximates  the  LHS  of  Eq.  (2.2).  To  expand  the  RHS  of  Eq.  (2.2)  in  terms  of  the  small  computational 
parameter 


(2.10) 


it  is  first  noted  that  fi*\xb,  t)  in  Eq.  (2.4)  can  be  expressed  as 

(2.11)  t)  =  fieq)(xh,  t)  +  wap^ea  ■  (ubf  -  uf) , 


so  that  the  RHS  of  Eq.  (2.2),  or  Eq.  (2.3),  can  be  rewritten  as 

fa{xf,  t)  «  fieq)(xf,  t)  +  (1  -  x)(l  -  1  /T)f£Hxf,  t)  +  wap^ea  •  ( XUbf-XUf  -  2 uw) 

3  3 

(2.12)  =  fieq){xf,  t)-{  1  -  x)(t  -  1  )Stwap-^ea  ■  V urea  +  wap-^ea  •  ( x«i>/  -  X«/  ~  2«u,) . 

Based  on  linear  interpolation,  Ubf  «  [(A  -  1  )uf  +  uw\/ A,  expanding  the  velocity  uw  at  Xf  near  the  wall 
(ecu,)  using  Taylor  series,  i.e.,  uw  =  Uf  +  A ^ea-Vu/,  one  obtains  Ubf  -  Uf  «  Stea^u.  Noticing  that 
(xb-Xf)  =  ea$t ,  and  equating  Eqs.  (2.9)  and  (2.12)  and  matching  terms  linear  in  5t  results  in  *  =  (2A-l)/r. 
For  Ubf  =  Uf ,  we  have  Ubf  —  Uf  =  0  in  Eq.  (2.12).  Matching  to  0(6t)  then  requires  x  —  (2A  -  1  )/(r  -  1). 
FH  found  that  Ubf  =  [(A  -  l)u/  +  uw\/ A  gives  computationally  stable  results  only  for  A  >  1/2.  Hence, 
they  proposed  that 

1  1  (2A  - 1)  1 

(2.13)  ubf  =  ^(A  -  1  )uf  +  —uw  ,  and  x  = - - - ,  for  A  >  -  , 


and 

(2A  -  1)  1 

(2.14)  ubf  —  uf ,  and  x  =  [y  >  for  A<2' 

To  recapitulate,  there  are  three  independent  assumptions  that  have  been  made  in  the  foregoing  deriva¬ 
tion.  They  are:  (i)  the  Chapman-Enskog  expansion  in  the  form  given  by  Eqs.  (1.9)  and  (1.10)  is  valid;  (ii) 
the  intrinsic  time  scale  of  the  unsteady  flow  must  be  much  large  compared  to  the  advection  time  on  the 
lattice  scale  given  by  Eq.  (2.5);  and  (iii)  the  lattice  spacing  must  be  small  compared  to  the  characteristic 
length  scale  of  the  flow  as  given  by  Eq.  (2.10)  so  that  the  Taylor  series  expansion  for  the  velocity  field  near 
the  wall  is  valid.  There  have  been  a  large  number  of  papers  in  the  existing  literature  regarding  the  validity 
and  usefulness  of  the  Chapman-Enskog  analysis  for  the  solution  of  the  Boltzmann  equation.  The  “slow  flow” 
condition  is  introduced  to  simplify  the  derivation  of  the  boundary  condition  for  fa ;  the  implication  of  this 
assumption  will  be  briefly  addressed  later  in  comparing  the  computational  results  with  that  based  on  the 
conventional  bounce-back  scheme.  The  last  assumption  is  a  typical  computational  resolution  requirement. 

Equation  (2.3)  is  essentially  a  linear  interpolation  (or  extrapolation)  and  is  used  continuously  in  the 
computation.  When  the  weighting  factor  x  becomes  too  large,  instability  may  develop.  For  1  >  A  >  1/2, 
|x|  =  |2A  —  l|/r  is  always  less  than  2  because  the  positivity  of  the  viscosity  in  the  LBE  scheme  requires 
r  >  1/2.  For  0  <  A  <  1/2,  |x|  =  |(2A  -  1  )/(r  -  1)|,  and  it  may  become  too  large  when  r  is  close  to  1.  To 
illustrate  this  point,  a  fully  developed  pressure  driven  2-D  channel  flow  is  considered.  The  grid  arrangement 
is  shown  in  Fig.  3.  For  steady  flow,  a  constant  pressure  gradient  Vp  along  the  z-direction  is  applied  and  can 
be  treated  as  a  body  force.  This  is  included  [23]  after  the  collision  step  by 


(2.15) 


fa(xi,  t)  =  fa{xi,  t)  +  wa5t^^ea-x 
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Fig.  3.  Lattice  configuration  in  channel  flow  simulations  with  an  arbitrary  A. 


where  x  is  the  unit  vector  along  the  x-axis.  The  boundary  condition  for  /a(xi,  t)  on  the  wall  follows  those 
given  by  Eqs.  (2.3),  (2.4),  (2.13),  and  (2.14).  At  the  inlet  ( i  =  1)  and  exit  (i  =  Nx ,  in  which  Nx  is  the 
number  of  lattices  in  the  x-direction)  the  following  is  imposed, 

(2.16)  fai}  —  1?  j)  —  fa( }  —  2,  j)  , 

(2.17)  =  NXi  j)  =  fa(i  ~  Nx  “  1)  j)  ♦ 

With  Eq.  (2.16),  the  velocity  profile  at  the  inlet,  ux(i  =  2,  j ),  is  not  needed.  Instead,  the  fully  developed 
velocity  profile  is  sought  as  part  of  the  solutions.  In  this  part  of  the  investigation,  Ny  =  35  is  used.  The  exact 
solution  for  the  velocity  profile  [given  by  Eq.  (3.1)]  is  used  as  the  velocity  initial  condition  which  differs  from 
the  final  steady  state  solution  due  to  numerical  errors.  The  equilibrium  distribution  function  fa  based  on 
the  exact  solution  for  the  velocity  profile  is  used  as  the  initial  condition  for  fa .  The  pressure  gradient  is  set 
to  ^  =  -1.0  x  10”6.  All  computations  are  carried  out  in  double  precision. 

For  A  <  1/2,  it  is  found  that  the  computation  is  unstable  for  certain  range  of  values  of  r.  In  Figure  4, 
solid  curves  are  the  stability-instability  boundaries  for  the  fully  developed  channel  flow  in  the  (A,  r)  space 
obtained  from  a  large  number  of  computations.  For  A  <  0.2,  the  computation  becomes  unstable  when  r  <  1. 
The  large  instability  region  is  an  apparent  source  of  concern  for  FH’s  scheme  when  A  <  1/2  since  lower 
viscosity  can  only  be  achieved  when  r  is  close  to  1/2. 

One  may  speculate  that  the  instability  in  the  above  example  results  from  the  lack  of  specifying  an  inlet 
velocity  profile,  ux(y ),  or  due  to  the  extrapolation  of  fa  at  the  inlet  given  by  Eq.  (2.16).  To  examine  this  pos¬ 
sibility,  a  channel  flow  entrance  problem  is  considered.  Uniform  velocity  profiles,  ux(y)  —  -  (H  / 12  pu)  (dp /dx) 
and  uy(y)  =  0  in  which  H  is  the  channel  height,  are  specified  at  i  —  1.5  (half  way  between  the  first  and 
second  lattices)  and  the  distribution  functions  fa(i  =  1,  j)  for  a  =  1,2,  and  8  are  obtained  using  Eq.  (2.3) 
with  x  =  0  in  accordance  with  A  =  1/2  at  i  =  1.5.  The  boundary  conditions  on  the  wall  are  based  on 
Eqs.  (2.3),  (2.4),  (2.13),  and  (2.14).  The  exit  boundary  condition  for  /a’s  is  given  by  Eq.  (2.17).  Hence  the 
extrapolation  for  fa  at  the  inlet  is  completely  eliminated  and  the  velocity  profiles  at  the  inlet  are  exactly 
given.  Two  types  of  initial  conditions  axe  used.  Whenever  possible,  the  equilibrium  distribution  functions 
corresponding  to  the  uniform  inlet  velocity  are  specified  at  t  =  0  throughout  the  flow  field.  This  works  for 
relatively  larger  values  of  r.  However,  instability  can  be  encountered  when  r  is  considerably  larger  than  the 
upper  solid  curve  shown  in  Fig.  4  for  the  same  value  of  A  (<  1/2).  A  second  type  of  initial  condition  is 
thus  implemented.  A  converged  solution  at  a  relatively  large  value  of  r  is  used  as  the  initial  condition  for  a 
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Fig.  4.  Regions  of  instability  in  the  LBE  computation  for  channel  flows  using  FH’s  boundary  condition ,  Eqs.  (2.3),  (2.4), 
(2.13),  and  (2.14),  with  A  <  1/2.  The  instability  region  for  the  fully  developed  2-D  channel  flow  (enclosed  by  solid  curves)  is 
slightly  smaller  than  the  instability  region  for  the  2-D  channel  flow  entrance  problem  ( enclosed  by  dashed  curves ). 


smaller  value  of  r.  The  value  of  r  is  incrementally  decreased  to  obtain  the  converged  solutions  for  the  new, 
smaller  values  of  r.  When  the  actual  instability  region  is  approached,  the  increment  in  r  is  maintained  as 
small  as  0.01  or  0.005.  In  the  computation,  &  =  —1.0  x  10“6,  Ny  =  35,  and  Nx  =  65  are  used.  When 
the  Reynolds  number  is  low  (due  to  the  use  of  the  small  pressure  gradient  and  larger  r),  the  exit  velocity 
profile  is  very  close  to  the  exact  solution  corresponding  to  the  fully  developed  channel  flow  which  validates 
the  solution  procedure. 

The  stability-instability  boundaries  for  the  channel  flow  entrance  problem  obtained  through  a  large 
number  of  computations  (dashed  curves)  are  also  shown  in  Fig.  4.  It  is  noted  that  the  stability-instability 
boundaries  for  the  channel  flow  entrance  problem  are  very  similar  to,  and  slightly  larger  than,  that  for  the 
fully  developed  channel  flow  despite  the  significant  difference  in  the  inlet  boundary  condition.  Thus  the 
source  of  the  instability  must  be  a  result  mainly  due  to  the  implementation  of  the  solid  wall  condition.  An 
alternative  scheme  must  be  developed  to  overcome  this  shortcoming. 

2.2.  Improved  treatment  for  curved  boundary.  We  realize  that  the  flexibility  in  the  construction 
of  is  the  key  to  achieve  an  improved  computational  stability  as  well  as  accuracy.  Since  x  =  (2  A  —  l)/(r  —  1) 
given  by  Eq.  (2.14)  leads  to  a  larger  value  of  x  when  r  is  close  to  1,  it  is  desirable  to  reduce  the  magnitude  of 
X  by  increasing  the  magnitude  of  the  denominator  in  the  expression  for  x*  For  A  >  1/2,  Ubf  is  the  fictitious 
fluid  velocity  inside  the  solid  and  the  denominator  for  x  is  r.  For  A  <  1/2,  Ubf  was  chosen  by  FH  to  be  Uf 
which  is  the  fluid  velocity  at  Xf  and  it  leads  to  (r  -  1)  in  the  denominator  for  x-  Thus,  we  propose  to  use 
Eq.  (2.13)  for  A  >  1/2  and  use 

(2.18)  Ubf  =  Uff  =  u(xf  4-  e&Sty  t) ,  for  A<i. 

Thus 


Ubf  -  Uf  =  u(xf  +  eaSt ,  t)  -  u(xf ,  t)  =  — JtV  Uf  -eoc . 


-t(1  -  x)(l  -  1/t)  -  X  =  2A  -  t 


(2.19) 

This  requires 

(2.20) 
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Fig.  5.  Regions  of  stability  and  instability  in  the  LBE  computation  for  channel  flow  using  the  present  boundary  condition, 
Eqs.  (2.3),  (2.4),  (2.13),  (2.18),  and  (2.21),  with  A  <  1/2,  similar  to  Fig.  4.  The  instability  region  for  the  fully  developed  2-D 
channel  flow  ( above  the  solid  curve)  is  almost  the  same  as  the  instability  region  for  the  2-D  channel  flow  entrance  problem 
( above  the  dashed  curve). 


to  match  the  0(St)  terms  in  equating  Eqs.  (2.9)  and  (2.12).  Hence 


(2.21) 


_  (2A-1) 

(r-2)  ’ 


for 


To  test  the  improvement  in  the  stability,  the  steady  state,  fully  developed,  pressure  driven  2-D  channel  flow 
is  again  considered.  Eqs.  (2.18)  and  (2.21)  are  used  in  lieu  of  Eq.  (2.14).  The  rest  of  the  implementation  is 
exactly  the  same  as  described  in  the  previous  section.  The  solid  curve  in  Fig.  5  shows  the  stability-instability 
boundary  in  the  (A,  r)  space  for  the  fully  developed  channel  flow.  By  comparing  Fig.  5  with  Fig.  4,  the 
improvement  in  the  stability  of  the  present  treatment  for  this  case  is  clearly  seen.  The  present  treatment 
moves  the  instability  region  around  r  =  1  caused  by  FH’s  boundary  condition  upward  in  the  parameter 
space  (A,  r)  to  the  region  around  r  =  2.  This  would  enable  us  to  use  r  in  the  interval  (0.5, 1.2]  for  small 
viscosity. 

For  the  channel  flow  entrance  problem,  boundary  conditions  at  the  inlet  and  exit  and  the  procedure 
for  specifying  the  initial  conditions  are  the  same  as  described  in  the  last  section.  Eqs.  (2.18)  and  (2.21) 
are  used  to  replace  Eq.  (2.14)  for  the  solid  wall.  The  stability-instability  boundary  in  the  (A,  r)  space  for 
the  entrance  flow  problem  is  shown  by  the  dashed  curve  in  Fig.  5.  Close  agreement  between  the  stability- 
instability  boundaries  for  the  two  cases  in  Fig.  5  suggests  that  the  improvement  in  the  computational  stability 
is  not  related  to  the  treatment  of  the  inlet  boundary  conditions.  The  improvement  results  rather  from  the 
different  treatment  in  the  solid  wall  boundary  condition.  One  advantage  of  the  present  treatment  of  solid 
wall  is  its  insensitivity  to  other  boundary  conditions  such  as  inlet  and  outlet  boundary  condition  in  channel 
flows,  as  indicated  by  the  results  shown  in  Fig.  5.  The  instability  in  the  present  treatment  is  mainly  due  the 
denominator  (r  —  2)  in  Eq.  (2.21).  A  direct  consequence  of  this  improvement  is  that  lower  values  of  r,  or 
lower  viscosity  v,  can  now  be  used. 

As  one  may  speculate  at  this  point  that  u(xf  4-  2 t)  can  also  be  used  for  u&/  when  A  <  1/2.  This 
would  further  improve  the  stability  since  x  =  (2  A  -  l)/(r  -  3).  This  is  correct  in  principle.  However,  since 
the  use  of  Uj  as  already  allows  the  use  of  r  whose  value  is  close  to  1/2,  there  is  little  practical  need  to 
use  Uf  that  is  too  far  away  from  the  wall. 
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For  transient  flows,  a  second-order  extrapolation  can  be  used  for 


-Uf  +  -£UW  +  A(1  +  A)  K 


(1  +  A)t if  +  A  u/f] 


(1-A). 


—  (1  +  A) A  'Uf+  A(H-A)Uu'’  - -  2' 

This  treatment  helps  to  improve  the  accuracy  in  the  velocity  approximation  when  u(x,  t)  is  not  well  resolved 
near  the  wall.  Finally,  it  is  easily  seen  that  the  present  boundary  condition  treatment  can  be  extended  to 
3-D  flow  problems  involving  curved  geometry.  The  efficacy  of  such  an  extension  will  be  examined  elsewhere. 

3.  Results  and  Discussions.  For  the  proposed  boundary  condition  treatment  to  be  useful,  several 
issues  need  to  be  addressed:  spatial  and  temporal  accuracy,  ability  to  handle  geometric  singularity,  and  the 
flexibility  to  handle  complex  geometry.  Channel  flows  with  constant  and  sinusoidally  oscillating  pressure 
gradients  with  analytic  solutions  are  used  to  assess  the  spatial  and  temporal  accuracy.  The  Stokes  first 
problem  (i.e.,  the  flow  due  to  an  impulsively  started  wall)  allows  one  to  examine  the  response  of  the  computed 
flow  field  to  an  imposed  singular  acceleration.  The  standard  lid-driven  cavity  flow  has  a  bounded  domain  but 
possesses  stress  or  vorticity  singularities  near  the  corners  between  the  moving  and  stationary  walls.  Finally, 
flow  over  a  column  of  circular  cylinders  is  the  case  used  to  assess  the  impact  of  the  boundary  treatment  on 
the  accuracy  of  the  flow  field  around  a  curved  boundary. 

3.1.  Pressure  driven  channel  flows.  At  steady  flow,  the  exact  solution  for  the  z- velocity  profile  is 
given  by 

.X  .  .  1  dpH2  ,  n  ^  ^  1 


for  A  > 


0<7?<  1, 


where  H  =  (Ny  -  3  +  2 A)  and  77  =  y/H  =  (j  -  2  +  A)/# .  To  assess  the  computational  error  of  the  LBE 
solution  of  the  velocity,  u(y),  the  following  relative  L2-norm  error  is  defined 


rH 

/  [u(y) -u0(y)]2dy 

Jo 


uo(y)dy 


With  the  oscillating  pressure  gradient,  ^  =  Be  lut ,  the  exact  solution  can  be  easily  expressed  in 
complex  variables.  An  important  parameter  in  this  flow  is  the  Stokes  number  St  defined  as, 


[6.6)  or  =  . 

y/v/UJ 

The  Stokes  number  is  the  ratio  of  the  channel  height  H  to  the  thickness  of  the  Stokes  layer  ^/vjw.  Since 
the  error  can  vary  with  time,  a  time  average  over  one  period  (T  =  27t/cj)  is  needed  and  the  relative  error  is 


[u(y,  t)-uo(y,  t)]2  dydt 


nH 

«o(2/>  t)dydt 

_ 


In  the  lattice  BGK  model,  6t  =  Sx  =  Sy  =  1.  Comparing  with  the  channel  height  H  =  (Ny  —  3  4-  2A),  the 
dimensionless  grid  size  (or  grid  resolution)  is  H. 


ll 
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Fig.  6.  Convergence  of  the  2-D  steady  state  pressure-driven  channel  flow  simulations  using  the  present  boundary  condition. 
The  symbols  represent  numerical  results  and  the  straight  lines  are  the  least-square  fit  of  the  data,  (a)  Convergence  of  global 
relative  L,2~norm  error  in  u(y)  on  H  =  ( Ny  —  3  +  2A).  (b)  Quadratic  convergence  of  the  slip  velocity  at  wall  uw.  (c)  Relative 
L2~norm  error  in  u{y)  as  a  function  of  A. 


Figure  6(a)  shows  the  dependence  of  the  relative  L2-norm  error  on  the  channel  height  H  for  r  =  0.55  and 
A  =  0.0,  0.25  and  0.5.  A  maximum  value  of  Ny  —  131  is  used.  The  second-order  accuracy  is  demonstrated 
in  the  range  of  H  investigated.  It  has  been  well  established  that  the  accuracy  of  the  LBE  method  for  the 
interior  points  is  of  second  order.  The  fact  that  the  overall  accuracy  is  of  second  order  in  the  present  case 
means  that  the  accuracy  in  the  boundary  condition  is  at  least  of  second  order.  It  is  worth  to  note  that  the 
derivation  given  in  Section  2  ensures  that  fa  is  second  order  accurate  near  the  wall.  It  does  not  guarantee 
the  second  order  accuracy  of  the  velocity  field  near  the  wall.  To  address  this  issue,  the  wall  slip  velocity, 
uw  =  ux(y  =  0),  is  evaluated  using  a  second  order  extrapolation  based  on  ux(y  —  A),  ux(y  =  1  +  A)  and 
ux(y  =  2  +  A).  Since  the  true  wall  velocity  in  the  pressure  driven  channel  flow  is  zero,  the  wall  slip  velocity 
uw  provides  a  measure  of  the  accuracy  for  the  treatment  of  the  wall  velocity.  Fig.  6(b)  shows  the  dependence 
of  uw ,  normalized  by  the  centerline  velocity  umax  =  -{H2 /8pv)(dp/dx),  on  H  for  A  =  0.0,  0.25,  and  0.5 
with  r  =  0.55.  Quadratic  convergence  is  clearly  observed  in  all  three  cases  which  demonstrate  the  second 
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order  accuracy  of  the  velocity  field  near  the  solid  wall.  This  is  entirely  consistent  with  the  results  shown 
in  Fig.  6(a)  which  involves  global  convergence  rather  than  the  local  ( y  =  0)  convergence.  Fig.  6(c)  shows 
the  relative  error  as  a  function  of  A  using  the  present  boundary  treatment  [Eqs.  (2.3),  (2.4),  (2.13),  (2.18), 
and  (2.21)]  for  0  <  A  <  1.  The  error  in  the  range  of  0  <  A  <  1/2  is  comparable  to  those  in  the  range 
of  1/2  <  A  <  1.  The  present  boundary  condition  treatment  does  not  induce  larger  computational  error 
and  is  substantially  more  robust.  Furthermore  the  2nd-order  accuracy  is  achieved  in  general  by  the  present 
treatment  for  A  <  1/2. 
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Fig.  7.  Dependence  of  the  1,2-norm  error  in  u(y)/um&x  on  the  lattice  resolution  H  =  {Ny  —  3  +  2A)  in  oscillating  pressure 
driven  channel  flow.  Stokes  number  St  =  ff/ yj vjw  =  1  and  8.  The  symbols  represent  numerical  results  and  the  straight  lines 
are  the  least-square  fit  of  the  data. 

Figure  7  shows  the  dependence  of  the  relative  I/2-norm  error  on  the  channel  height  H  in  the  oscillating 
pressure  driven  channel  flow  for  Stokes  number  St  =  —  1  and  8.  For  St  =  1,  the  Stokes  layer  is 

as  thick  as  the  channel  height  H.  For  A  =  0.25,  0.5,  and  0.75,  second-order  accuracy  in  space  is  clearly 
demonstrated.  Since  the  time  step  St  in  LBE  is  equal  to  the  spatial  resolution  8Xi  the  accuracy  in  time  must 
also  be  of  second-order  in  order  for  the  time-averaged  L2-norm  error  to  have  a  slope  of  2.0  in  Fig.  7.  For 
St  =  8,  the  Stokes  layer  thickness  is  about  1/8  of  the  channel  height  so  that  the  computational  error  due 
to  the  insufficient  resolution  of  the  Stokes  layer  is  a  significant  part  of  the  error.  For  A  =  0.25,  the  first 
lattice  in  the  flow  field  is  only  a  quarter  of  the  lattice  size  away  from  the  wall.  The  Stokes  layer  is  thus 
better  resolved  for  A  =  0.25  (denoted  by  solid  circles  in  Fig.  7)  than  for  A  =  0.5  and  0.75.  However,  as  H 
increases,  the  difference  between,  A  =  0.25  and  A  =  0.5  and  0.75  becomes  smaller  since  all  have  reasonable 
resolutions  in  the  Stokes  layer.  Although  the  slope  for  the  error  curve  for  A  =  0.25  is  observed  to  be  about 
1.5  that  is  less  than  2,  it  is  an  indication  of  the  better-than-expected  accuracy  at  the  low  resolution  end. 

3.2.  Stokes  first  problem:  Flow  due  to  an  impulsively  started  wall.  For  a  wall  located  at  y  =  0 
that  is  impulsively  started,  an  unsteady  Stokes  layer  of  thickness  0(y/vi)  develops  near  the  wall.  For  a 
fixed-grid  computation,  the  error  at  small  time  is  expected  to  be  large  due  to  insufficient  spatial  resolution. 
In  the  LBE  method,  this  is  also  compounded  by  the  use  of  fixed  6t  (=  Sx  —  Sy  =  1).  Figure  8  shows  the 
velocity  profiles  at  t  =  100  (in  lattice  unit).  The  wall  velocity  is  U  =  0.1  in  lattice  unit.  The  relaxation 
time  r  =  0.52  gives  kinematic  viscosity  v  —  0.0067.  Similar  to  the  oscillating  pressure  driven  channel  flow, 
the  error  is  smaller  for  A  =  0.25  than  for  A  =  0.5  and  0.75  due  to  a  better  spatial  resolution  near  the  wall. 
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Fig.  8.  Velocity  profiles  u(y ,  t)/U  of  an  impulsively  started  plate  (Stocks  first  problem)  with  A  —  0.25,  0.5,  and  0.75 
(symbols),  att—  100  (in  lattice  unit).  The  solid  line  represents  the  exact  solution  of  the  problem.  The  bounce  back  on  the  link 
(BBL)  always  sets  A  =  0.5. 


Fig.  9.  Relative  L2-norm  error  of  the  velocity  profile  u(y,  t)  during  the  initial  transient  of  the  impulsively  started  plate 
with  various  values  of  A  (0.25,  0.5,  and  0.75).  The  “ linear ”  version  of  the  boundary  condition  corresponds  to  Eq.  (2.13).  The 
“ quadratic ”  version  corresponds  to  Eq.  (2.22).  BBL  is  limited  to  A  =  0.5  only. 


Figure  9  shows  the  temporal  variation  of  the  relative  I/2-norm  error  defined  as 

V2 


(3.5) 


E2(t)  = 


|  J  [«(y,  t)  -  u0(y,  i)]2  dy  j 

"  r 


Uo(y,  t)dy 


for  A  =  0.25,  0.5,  and  0.75.  The  result  using  the  standard  bounce-back  on  the  link  (BBL)  scheme,  which 
always  sets  A  =  0.5,  is  also  shown.  The  large  relative  errors  in  the  beginning  are  due  to  the  smaller  values 
of  the  denominator  in  the  above  equation.  It  should  be  emphasized  that  this  flow  at  small  time  is  difficult 
to  deal  with  for  any  computational  technique  due  to  the  singular  acceleration  and  large  spatial  gradient. 
For  an  impulsively  started  Couette  flow,  the  long-time  solution  approaches  the  exact  linear  velocity  profile 
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because  the  LBE  method  is  a  second-order  accurate  one.  It  is  interesting  to  note  that  the  present  boundary 
condition  treatment  for  A  =  0.5  gives  a  slightly  smaller  error  than  the  BBL  scheme  in  this  highly  transient 
case.  In  such  a  transient  flow,  the  computational  accuracy  in  the  near-wall  region  is  typically  dictated  by 
the  near-wall  spatial  resolution  which  must  be  smaller  than  the  Stokes  layer  thickness  in  order  to  resolve 
the  local  flow  field.  In  a  finite  difference  calculation  for  such  a  flow,  St  and  Sx  can  be  independently  chosen. 
If  Sx  is  not  sufficiently  small,  further  reduction  in  St  will  not  lead  to  improvement  in  accuracy.  At  small  £, 
neither  the  BBL  scheme  nor  the  present  treatment  resolved  the  Stokes  layer  so  that  the  error  is  large.  After 
the  Stokes  layer  grows  to  certain  thickness,  the  spatial  resolution  becomes  adequate  and  the  accuracy  then 
improves.  In  view  of  the  “slow  flow”  condition  (2.5)  introduced  in  the  derivation,  the  performance  of  the 
current  boundary  treatment  is  comparable  or  better  than  the  conventional  bounce-back  on  the  link  scheme. 


Fig.  10.  Velocity  profiles  at  the  center  (x/H  —  1/2)  in  lid-driven  cavity  flow  with  various  values  of  A  (0.1,  0.5,  and  0.9) 
at  Re  —  100. 


3.3.  Flow  in  a  lid-driven  square  cavity.  Figure  10  shows  the  velocity  profiles  at  the  center  ( x/H  = 
1/2)  of  the  cavity  of  width  H  at  Re  =  100  with  r  =  0.6.  Only  35  x  35  lattices  are  used  and  the  cavity  width 
is  H  =  ( Nx  —  3  +  2 A)  =  32  -f  2 A.  This  requires  the  lid  velocity  to  be  U  =  vRe/H  =  3.33 /H  in  the  lattice 
unit.  It  has  a  negligible  compressibility  effect  for  H  ~  32.  A  well-resolved  finite  difference  solution  for  the 
velocity  field  based  on  the  stream  function- vorticity  formulation  is  also  shown  for  comparison.  The  velocity 
profile  with  A  =  0.1  agrees  well  the  finite  difference  solution.  For  A  =  0.5,  the  result  is  rather  reasonable 
with  such  a  resolution.  The  difference  is  slightly  larger  on  the  negative  velocity  part  for  A  =  0.9.  The 
corner  singularity  in  stress  (or  vorticity)  is  well  handled  for  r  =  0.6  and  Nx  =  35.  However  for  r  close  to  0.5 
and  with  Nx  =  35,  the  corner  singularity  induces  wiggles  in  the  velocity  field.  This  issue  will  be  examined 
elsewhere.  Flow  field  for  Re  =  1000  is  obtained  with  67  x  67  lattices  using  A  =  0.1,  0.5,  and  0.9.  Similar 
behavior  in  the  velocity  profiles  is  observed. 

3.4.  Uniform  flow  over  a  column  of  circular  cylinders.  To  simulate  the  external  flow  over  a 
single  cylinder  would  require  placing  the  outer  boundary  far  away  from  the  cylinder.  In  order  to  keep  the 
computational  effort  at  a  reasonable  level  in  using  constant  space  lattices,  a  column  of  circular  cylinders  of 
radius  r  and  center-to-center  distance  H  is  considered  instead.  The  flow  field  that  needs  to  be  computed  is 
thus  limited  to  —H  <  y  <  H.  At  y  =  —H,  the  lattice  is  j  =  2.  The  boundary  conditions  at  j  =  1  for  fa  s 
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Fig.  11.  Flow  past  a  column  of  cylinders.  Velocity  profiles  at  x  —  0  for  uniform  flow  over  a  column  of  cylinders.  The 
cylinder  has  a  diameter  (2 r)  of  7  lattice  units.  The  cylinder  center-to-center  distance  H  =  70  in  lattice  units. 

are  given  by  the  following  symmetry  properties, 

fo(i,  1)  =  Mh  3),  Mh  1)  =  fi(h  3),  h(i ,  1)  =  Mh  3), 

(3.6)  Mh  1)  =  Mh  3),  Mh  1)  =  Mh  3),  Mh  1)  =  Mh  3), 

Mh  1)  =  Mh  3),  Mh  1)  =  Mh  3),  Mh  1)  =  Mh  3). 

Similar  conditions  hold  at  y  =  H  for  j  =  Ny.  At  the  inlet,  the  uniform  velocity,  u(y)  =  U,  is  specified  at 
i  =  1.5.  Using  A  =  0.5,  x  —  0,  Eq.  (2.3)  is  applied  to  obtain  fa' s  at  i  =  1.  At  the  exit,  a  simple  extrapolation 
is  used, 

(3.7)  fa(Nx ,  j)  =  2fa(Nx  -  1,  j)  -  fa(Nx  -  2,  j) ,  for  a  =  4,  5,  and  6  . 

On  the  surface  of  the  circular  cylinder,  Eqs.  (2.3),  (2.4),  (2.13),  (2.18),  and  (2.21)  proposed  in  this  paper  are 
used  to  update  the  boundary  conditions  for  /a’s. 

Figure  11  shows  the  velocity  profile  u(x  —  0,  y)/U  for  H/r  =  20  at  Re  =  2 Ur/u  =  10  using  r  =  3.5. 
Two  values  of  relaxation  time  r,  0-505  and  0.525,  are  used.  For  r  =  3.5,  there  are  only  7  lattices  from  the 
front  to  the  back  stagnation  points.  The  finite  difference  solution  is  obtained  using  u-ip  formulation  with 
body-fitted  coordinates  [26]  and  over  200  grid  points  are  distributed  along  the  upper  half  of  the  circle.  These 
two  solution  with  r  =  0.505  and  r  =  0.525  are  virtually  identical  to  each  other  and  they  are  both  close  to  the 
finite  difference  solution.  Fig.  12  shows  the  centerline  (y  =  0)  velocity  variations,  upstream  and  downstream, 
at  Re  =  10  and  40.  The  sharp  gradient  near  the  front  stagnation  point,  the  length  of  the  separation  bubble, 
the  maximum  of  the  separation  bubble  velocity,  and  the  recovery  of  the  wake  velocity  are  all  in  excellent 
agreement  with  the  well  resolved  finite  difference  solution. 

As  can  be  seen  now  that  an  important  improvement  of  the  present  boundary  condition  treatment  over 
the  bounce-back  scheme  is  that  it  can  preserve  the  accuracy  of  the  geometry  under  consideration.  To  further 
demonstrate  this  point,  consider  flow  over  a  circular  cylinder  of  radius  r  with  the  coordinate  centered  at  the 
center  of  the  cylinder.  For  r  =  3.4  and  3.8,  the  front  stagnation  points  are  located  at  x  =  -3.4  and  -3.8, 
respectively.  With  the  bounce-back  on  the  link  (BBL)  scheme,  the  front  stagnation  points  in  both  cases  will 
be  placed  at  x  =  -3.5  which  is  half-way  between  the  lattice  at  x  =  -4  and  x  =  -3  on  the  centerline.  In  the 
present  method,  A  =  0.6  and  0.2  for  r  —  3.4  and  3.8,  respectively.  The  difference  in  A  can  be  accurately 
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Fig.  12.  Flow  past  a  column  of  cylinders.  Centerline  ( y  =  0)  velocity  variation  for  a  uniform  flow  over  a  column  of 
cylinders  for  Re  =  10  and  40.  The  cylinder  diameter  2 r  is  only  of  7  lattice  units.  Finite  difference  results  are  based  on  u-ip 
formulation  and  are  well  resolved. 


Fig.  13.  Flow  past  a  2-D  cylinder.  Similar  to  Fig.  11.  Comparison  of  the  velocity  profiles  at  x  =  0  for  the  radius  r  =  3.0, 
3.2,  3.4,  3.5,  3.6,  3.8  and  4.0  with  Re  =  10,  r  =  0.7,  and  H/r  =  20. 


incorporated  in  the  evaluation  of  /5(x& ,  t).  This  implies  that  although  the  boundary  links  for  r  =  3.4  will  be 
different  from  those  for  r  =  3.8,  the  flow  fields  based  on  r  =  3.4  and  r  =  3.8  should  be  nearly  the  same  when 
the  coordinates  are  normalized  by  the  radius  r .  To  validate  this  point,  a  series  of  computations  are  carried 
out  for  r  —  3.0,  3.2,  3.4,  3.5,  3.6,  3.8  and  4.0  for  H/r  —  20  at  Re  =  10.  The  profiles  of  the  dimensionless 
rr-component  velocity  ux(x ,  y)/U  as  a  function  of  y/r  at  x  =  0  are  compared  for  these  seven  different  radii 
r  in  Fig.  13.  Excellent  agreement  is  observed.  Figure  14  compares  the  ux(x ,  y)/U  as  a  function  of  x/r  at 
y  —  0  for  both  the  downstream  and  upstream  regions  for  these  seven  different  radii.  Again,  all  seven  cases 
compare  very  well  even  in  the  near  wall  region.  Furthermore,  the  LBE  solutions  with  different  radii  and 
the  same  Reynolds  number  all  well  agree  with  the  u )  —  <j)  finite  difference  solution  of  the  incompressible  NS 
equation  with  a  well  resolved  resolution.  This  clearly  demonstrates  that  the  present  boundary  condition 
treatment  has  maintained  geometric  fidelity  even  with  coarse  grid  resolutions. 
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Fig.  14.  Flow  past  a  2-D  cylinder.  Similar  to  Fig.  12.  Comparison  of  the  centerline  (y  —  0)  velocity  for  radius  r  =  3.0, 
3.2,  3.4,  3.5,  3.6,  3.8  and  4.0  with  Re  =  10,  r  =  0.7,  and  H/r  =  20. 


It  is  noted  that  the  interpolation  for  fa(&b,  t)  given  by  Eq.  (2.3)  is  carried  along  the  line  in  the  direction 
of  ea.  The  results  for  flow  over  a  cylinder  are  quite  satisfactory.  Other  interpolation  procedures  can  certainly 
be  devised  to  use  more  information  on  neighboring  lattices  in  the  flow  field.  However,  this  will  result  in  a  lot 
more  complications  in  the  implementation.  It  is  not  clear  if  such  an  attempt  will  necessarily  lead  to  further 
improvement  over  the  present  approach. 

4.  Conclusion.  In  this  work  a  second-order  accurate  boundary  condition  treatment  for  the  lattice 
Boltzmann  equation  is  proposed.  A  series  of  studies  are  conducted  to  systematically  validate  the  accuracy 
and  examine  the  robustness  of  the  proposed  boundary  condition  in  steady  and  unsteady  flows  involving 
flat  and  curved  walls.  Compared  with  the  existing  method  for  treating  boundary  condition  in  the  lattice 
Boltzmann  method,  the  proposed  treatment  has  the  following  advantages,  (i)  It  can  preserve  the  geometry 
of  interest  without  truncating  it  into  a  series  of  stair  steps,  (ii)  The  boundary  treatment  generally  results  in 
solutions  of  second-order  accuracy  for  the  velocity  field  in  space,  and  in  time  for  some  cases,  (iii)  Compared 
with  the  widely  used  bounce-back  on  the  link  scheme,  the  present  treatment  gives  comparable  or  better 
results  for  the  flow  field  under  otherwise  identical  computational  parameters. 
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